Preparations

Load the necessary libraries

library(car)
library(rstanarm) # for fitting models in STAN
library(brms) # for fitting models in STAN
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(rstan) # for interfacing with STAN
library(emmeans) # for marginal means etc
library(DHARMa) # for residual diagnostics
library(broom) # for tidying outputs
library(tidybayes) # for more tidying outputs
library(ggeffects) # for partial plots
library(broom.mixed) # for summarising models
library(ggeffects) # for partial effects plots
library(tidyverse) # for data wrangling etc
library(patchwork) # for multiple plots

Scenario

Loyn (1987) modelled the abundance of forest birds with six predictor variables (patch area, distance to nearest patch, distance to nearest larger patch, grazing intensity, altitude and years since the patch had been isolated).

Regent honeyeater

Format of loyn.csv data file

ABUND DIST LDIST AREA GRAZE ALT YR.ISOL
.. .. .. .. .. .. ..
ABUND Abundance of forest birds in patch- response variable
DIST Distance to nearest patch - predictor variable
LDIST Distance to nearest larger patch - predictor variable
AREA Size of the patch - predictor variable
GRAZE Grazing intensity (1 to 5, representing light to heavy) - predictor variable
ALT Altitude - predictor variable
YR.ISOL Number of years since the patch was isolated - predictor variable

The aim of the analysis is to investigate the effects of a range of predictors on the abundance of forest birds.

Read in the data

loyn <- read_csv("../public/data/loyn.csv", trim_ws = TRUE)
glimpse(loyn)
## Rows: 56
## Columns: 7
## $ ABUND   <dbl> 5.3, 2.0, 1.5, 17.1, 13.8, 14.1, 3.8, 2.2, 3.3, 3.0, 27.6, 1.8…
## $ AREA    <dbl> 0.1, 0.5, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.…
## $ YR.ISOL <dbl> 1968, 1920, 1900, 1966, 1918, 1965, 1955, 1920, 1965, 1900, 19…
## $ DIST    <dbl> 39, 234, 104, 66, 246, 234, 467, 284, 156, 311, 66, 93, 39, 40…
## $ LDIST   <dbl> 39, 234, 311, 66, 246, 285, 467, 1829, 156, 571, 332, 93, 39, …
## $ GRAZE   <dbl> 2, 5, 5, 3, 5, 3, 5, 5, 4, 5, 3, 5, 2, 1, 5, 5, 3, 3, 3, 2, 2,…
## $ ALT     <dbl> 160, 60, 140, 160, 140, 130, 90, 60, 130, 130, 210, 160, 210, …

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ log(\mu_i) = \boldsymbol{\beta} \bf{X_i} \]

where \(\boldsymbol{\beta}\) is a vector of effects parameters and \(\bf{X}\) is a model matrix representing the additive effects of the scaled versions of distance (ln), distance to the nearest large patch (ln), patch area (ln), grazing intensity, year of isolation and altitude on the abundance of forest birds.

Fit the model

MCMC sampling diagnostics

Model validation

Partial effects plots

Model investigation

Further analyses

Contrasts

brms

loyn.list <- with(loyn, list(AREA = c(min(AREA), mean(AREA), max(AREA))))

newdata <- loyn.brm4b %>%
  emmeans(~ fGRAZE | AREA, at = loyn.list, type = "response") %>%
  pairs() %>%
  as.data.frame()
head(newdata)
loyn.brm4b %>%
  emmeans(~ fGRAZE | AREA, at = loyn.list) %>%
  regrid() %>%
  pairs()
## AREA =    0.1:
##  contrast estimate lower.HPD upper.HPD
##  1 - 2       9.294     0.694     19.11
##  1 - 3      10.323     0.894     19.59
##  1 - 4      16.191     7.431     25.33
##  1 - 5      17.065     9.666     25.87
##  2 - 3       0.966    -6.671      9.07
##  2 - 4       6.870    -0.534     14.44
##  2 - 5       7.736     1.068     14.53
##  3 - 4       5.915    -1.843     12.67
##  3 - 5       6.682     1.288     13.77
##  4 - 5       0.722    -3.801      7.83
## 
## AREA =   69.3:
##  contrast estimate lower.HPD upper.HPD
##  1 - 2      -3.183   -10.440      4.18
##  1 - 3      -5.365   -13.947      4.42
##  1 - 4     -24.992   -71.614      1.07
##  1 - 5       8.318   -11.813     23.73
##  2 - 3      -2.187   -13.059      9.20
##  2 - 4     -21.617   -66.587      7.96
##  2 - 5      11.241   -10.116     29.02
##  3 - 4     -19.818   -64.696     12.25
##  3 - 5      13.350    -9.618     30.33
##  4 - 5      32.628    -9.247     78.21
## 
## AREA = 1771.0:
##  contrast estimate lower.HPD upper.HPD
##  1 - 2     -21.710   -57.156      5.77
##  1 - 3     -30.583   -75.732      9.54
##  1 - 4    -203.249 -1151.671     10.59
##  1 - 5     -31.248  -308.030     36.73
##  2 - 3      -9.001   -62.776     46.04
##  2 - 4    -178.997 -1144.401     46.68
##  2 - 5     -10.066  -286.823    102.55
##  3 - 4    -171.191 -1140.089     68.90
##  3 - 5      -0.977  -282.146    125.12
##  4 - 5     155.627  -356.586   1343.55
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95
newdata <- loyn.brm4b %>%
  emmeans(~ fGRAZE | AREA, at = loyn.list, type = "response") %>%
  pairs() %>%
  gather_emmeans_draws()

newdata %>%
  median_hdci() %>%
  ggplot() +
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_pointrange(aes(y = .value, ymin = .lower, ymax = .upper, x = contrast)) +
  facet_wrap(~AREA) +
  coord_flip()

emmeans(loyn.brm4b, ~ fGRAZE | AREA, at = loyn.list, type = "response") %>%
  gather_emmeans_draws()
newdata.p <- newdata %>% summarise(P = sum(.value > 1) / n())
g <- newdata %>%
  ggplot() +
  geom_vline(xintercept = 1, linetype = "dashed") +
  stat_slab(aes(
    x = .value, y = contrast,
    fill = stat(ggdist::cut_cdf_qi(cdf,
      .width = c(0.5, 0.8, 0.95),
      labels = scales::percent_format()
    ))
  ), color = "black") +
  scale_fill_brewer("Interval", direction = -1, na.translate = FALSE) +
  facet_grid(~ round(AREA, 1))

g + geom_text(data = newdata.p, aes(y = contrast, x = 1, label = round(P, 3)))

g + geom_text(data = newdata.p, aes(y = contrast, x = 1, label = paste("P = ", round(P, 3))), hjust = -0.2, position = position_nudge(y = 0.5))

Summary figure

References

Loyn, R. H. 1987. “Nature Conservation: The Role of Remnants of Native Vegetation.” In, edited by D. A. Saunders, G. W. Arnold, A. A. Burbridge, and A. J. M. Hopkins. Chipping Norton, NSW: Surrey Beatty & Sons.